//      sparse_vec.cc
//      
//      Copyright 2012 tqlong <tqlong@espada>
//      
//      This program is free software; you can redistribute it and/or modify
//      it under the terms of the GNU General Public License as published by
//      the Free Software Foundation; either version 2 of the License, or
//      (at your option) any later version.
//      
//      This program is distributed in the hope that it will be useful,
//      but WITHOUT ANY WARRANTY; without even the implied warranty of
//      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//      GNU General Public License for more details.
//      
//      You should have received a copy of the GNU General Public License
//      along with this program; if not, write to the Free Software
//      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
//      MA 02110-1301, USA.

#include "sparse_vec.h"

const SparseVector& SparseVector::operator += (const SparseVector& b)
{
  SparseVector a = *this;
  *this = a + b;
  return *this;
}

const SparseVector& SparseVector::operator -= (const SparseVector& b)
{
  SparseVector a = *this;
  *this = a - b;
  return *this;
}

std::ostream& operator << (std::ostream& os, const SparseVector& x)
{
  for (int i = 0; i < sz(x); i++) {
    os << "\t" << x.nz_[i] << ":" << x.at(i);
  }
  os << ENDL;
  return os;
}

SparseVector operator / (const SparseVector& x, double s)
{
  SparseVector tmp = x;
  for (int i = 0; i < sz(tmp); i++)
    tmp[i] /= s;
  return tmp;
}

SparseVector operator * (const SparseVector& x, double s)
{
  SparseVector tmp = x;
  for (int i = 0; i < sz(tmp); i++)
    tmp[i] *= s;
  return tmp;
}

SparseVector operator + (const SparseVector& x, const SparseVector& y)
{
  SparseVector tmp;
  int i = 0, j = 0;
  for (; i < sz(x) || j < sz(y); ) {
    if (j >= sz(y)) {
      tmp.pb(x[i]);
      tmp.nz_.pb(x.nz_[i]);
      i++;
    }
    else if (i >= sz(x)) {
      tmp.pb(y[j]);
      tmp.nz_.pb(y.nz_[j]);
      j++;
    }
    else if (x.nz_[i] < y.nz_[j]) {
      tmp.pb(x[i]);
      tmp.nz_.pb(x.nz_[i]);
      i++;
    }
    else if (x.nz_[i] > y.nz_[j]) {
      tmp.pb(y[j]);
      tmp.nz_.pb(y.nz_[j]);
      j++;
    }
    else {
      tmp.pb(x[i]+y[j]);
      tmp.nz_.pb(y.nz_[j]);
      i++; j++;
    }
  }
  return tmp;
}

SparseVector operator - (const SparseVector& x, const SparseVector& y)
{
  SparseVector tmp;
  int i = 0, j = 0;
  for (; i < sz(x) || j < sz(y); ) {
    if (j >= sz(y)) {
      tmp.pb(x[i]);
      tmp.nz_.pb(x.nz_[i]);
      i++;
    }
    else if (i >= sz(x)) {
      tmp.pb(-y[j]);
      tmp.nz_.pb(y.nz_[j]);
      j++;
    }
    else if (x.nz_[i] < y.nz_[j]) {
      tmp.pb(x[i]);
      tmp.nz_.pb(x.nz_[i]);
      i++;
    }
    else if (x.nz_[i] > y.nz_[j]) {
      tmp.pb(-y[j]);
      tmp.nz_.pb(y.nz_[j]);
      j++;
    }
    else {
      tmp.pb(x[i]-y[j]);
      tmp.nz_.pb(y.nz_[j]);
      i++; j++;
    }
  }
  return tmp;
}

double arma::norm(const SparseVector& x, int type)
{
  double s = 0;
  
  switch (type) {
  case 2: 
    for (int i = 0; i < sz(x); i++) s += sqr(x[i]);
    return ::sqrt(s);
  default:
    FATAL_ERROR ("norm " << type << " not supported");
  }
}

double arma::dot(const SparseVector& x, const SparseVector& y)
{
  double s = 0;
  int i = 0, j = 0;
  for (; i < sz(x) && j < sz(y); ) {
    if (x.nz_[i] < y.nz_[j]) {
      i++;
    }
    else if (x.nz_[i] > y.nz_[j]) {
      j++;
    }
    else {
      s += x[i] * y[j];
      i++; j++;
    }
  }
  return s;
}
